#include <cstdlib>
#include <cstdio>
#include <cmath>
#include "ising.h"

ising::ising (uint32 len, uint8 dim, double h, double j, double t)
{
	this->len = len;
	this->dim = dim;
	this->h = h;
	this->j = j;
	this->t = t;
	n_sites = int_pow (len, dim);

	n_parity = 2;

	lat = new uint8 [n_sites];
	// Randomize the lattice
	for (int i=0; i<(int)n_sites; i++)
		lat [i] = (rand () % 2) ? DOWN : UP;
}

ising::~ising ()
{
	delete lat;
}

void ising::print_lattice ()
{
	for (int i=0; i<(int)n_sites; ++i)
	{
		printf ("%d ", lat[i]);
		if ( !((i+1)%len) ) printf ("\n");
	}
}
void ising::flip_site (uint32 site) 
{
	if (site > n_sites) printf ("Site %u out of range.\n", site);
	else lat[site] = (lat[site] ? DOWN : UP);
}

int ising::spin_count (uint32 site, uint8 state)
{
	/*
	 * Count the total S_i \dot S_j for some site.
	 */
	static int i, dir;
	// Total S_i \dot S_j
	int spins = 0;
	int len_d = 1;

	for (i=0; i<dim; ++i)
	{
		for (dir=-1; dir<=1; dir+=2)
		{
			spins += 
				(
				(state ^
				lat [ (site+dir*len_d)%(n_sites) ]
				// if anti-aligned, return -1 = s\dot s_j
				) ? -1 : 1);
		}		
		// Shifting over to the next dimension/coordinate
		len_d *= len;
	}
	return spins;
}

double ising::get_energy (uint32 site, uint8 state)
{
	return -j*(spin_count (site, state)) - h*(state?1:-1);
}

void ising::random_update ()
{
	static int i, r;
	static double Ei, Ef;
	Ei = Ef = 0;

	for (i=0; i<(int)n_sites; ++i)
	{
		r = rand () % n_sites;
		Ei = get_energy (r, lat[r]);
		Ef = get_energy (r, lat[r] ? DOWN : UP);
		if (Ef < Ei)
			flip_site (r);
		else if (rand01() < exp(-(Ef-Ei)/t))
			flip_site (r);
	}
}

int ising::check_parity (int site)
{
	static int acc, len_d, i, coord;
	acc = 0; 
	len_d = 1;
	for (i=0; i<dim; ++i)
	{
		coord = (i / len_d) % len;
		acc += coord;
		len_d *= len;
	}
	return acc % n_parity;
}

void ising::checkerboard_update (int parity)
{
	static int i;
	static double Ei, Ef;
	Ei = Ef = 0;

	for (i=0; i<(int)n_sites; ++i)
	{
		if (check_parity (i) == parity)
		{
			Ei = get_energy (i, lat[i]);
			Ef = get_energy (i, lat[i] ? DOWN : UP);
			if (Ef < Ei)
				flip_site (i);
			else if (rand01() < exp(-(Ef-Ei)/t))
				flip_site (i);
		}
	}
}

void ising::swap_mag_field ()
{ 
	h = -h;
}

void ising::step (int method)
{
	switch (method)
	{
		case 0:
			for (int i=0; i<n_parity; ++i)
				checkerboard_update (i);
			break;
		case 1:
			random_update ();
			break;
		case 2:
			for (int i=0; i<n_parity; ++i)
				checkerboard_update ( rand() % n_parity );
		default:
			printf ("Invalid method chosen in step. Exiting.\n");
			exit (1);
	}
}
